Shock and awe with R
(Remi/David)
library(sf)
library(marmap)
library(raster)
library(tidyverse)
library(ggplot2)IMAGINE!
Let’s imagine that you are interested in a species in a given area and wish to know as much as possible about it. But, you can’t go out in the field because funding is running short. What you do have, however, is a certain knowledge of the tools that are available to you and you are wondering how far they can take you. Let’s find out.
Defining species and area of interest
# Let's select Atlantic cod as our species of interest
sp <- 'Gadus morhua'
# And the gulf of St. Lawrence as our study area
#define the corners
latmax <- 52.01312
latmin <- 45.52399
lonmax <- -55.73636
lonmin <- -71.06333
#
# latmax <- 46.25
# latmin <- 45.5
# lonmax <- -61.25
# lonmin <- -62.25Describe your species
Fishbase
We’ll start with a description of the species. First, let’s see what fishbase has to offer. This online data repository, along with sealifebase, contains a lot of information on marine and aquatic species all over the world and is accessible through the package [rfishbase](https://github.com/ropensci/rfishbase)
# Let's have a look at the information available our species ecology
ecol <- rfishbase::ecology(sp)
ecol <- cbind(colnames(ecol), t(ecol))
rownames(ecol) <- NULL
ecol <- ecol[ecol[,2] != 0, ] # remove 0
ecol <- ecol[!is.na(ecol[,2]), ] # remove NAs
knitr::kable(ecol, col.names = c('Descriptors', 'Attributes'))| Descriptors | Attributes |
|---|---|
| autoctr | 33 |
| sciname | Gadus morhua |
| StockCode | 79 |
| EcologyRefNo | 1371 |
| HabitatsRef | 1371 |
| Neritic | -1 |
| Intertidal | -1 |
| Oceanic | -1 |
| Estuaries | -1 |
| Herbivory2 | mainly animals (troph. 2.8 and up) |
| HerbivoryRef | 5743 |
| FeedingType | hunting macrofauna (predator) |
| FeedingTypeRef | 5743 |
| DietTroph | 4.09 |
| DietSeTroph | 0.179 |
| DietTLu | 4.34 |
| DietseTLu | 0.72 |
| DietRemark | Troph of adults from 7 studies. |
| DietRef | 26813 |
| FoodTroph | 4.29 |
| FoodSeTroph | 1 |
| FoodRemark | Trophic level estimated from a number of food items using a randomized resampling routine. |
| AddRems | Opportunistic predator that forages mainly at dawn and dusk (Refs. 1371, 46189). Larvae feed mainly on zooplankton while juveniles prey predominantly on benthic crustaceans; adults feed mainly on zoobenthos and fish (Refs. 5743, 9604, 26813) including juvenile cod. Fish prey becomes more common in the diet with increasing body size (Refs. 1371, 89387). Adults may cover large distances during the feeding period (Ref. 89387). |
| Young cod are also pre | yed upon by different fish species and octopus. Adult cod are prey items of top predators like sharks, rays, whales, dolphins, seals, and sea birds (Refs. 9023, 9581, 26954, 43651, 45735). |
| In the Baltic it grows | up to 5 kg weight in 7-8 years; in the North Sea it reaches 8 kg in the same time span . Natural mortality for adults of both stocks is assumed to be around M=0.2, resulting in a mean adult life expectancy and mean duration of the reproductive phase of 5 years (Ref. 88171). |
| Parasites of the speci | es include protozoans (trypanosome), myxosporidians, monogeneid, trematodes, cestodes, nematodes, acanthocephalan, hirudinid and copepods (Ref. 5951). |
| Schooling | -1 |
| SchoolingFrequency | sometimes |
| SchoolingLifestage | juveniles and adults |
| SchoolShoalRef | 1371 |
| AssociationsRemarks | Generally considered a demersal fish although its habitat may become pelagic under certain hydrogrphic conditions, when feeding or spawning. There is some evidence that cod leave the bottom and school pelagically to spawn in preferred temperatures when bottom tempetatures are unsuitable. Gregarious during the day, forming compact schools that swim between 30-80 m above the bottom, and scatter at night (Ref. 1371). Schooling behavior may be adaptive for feeding. Reproductive behavior during spawning involves the circling of a female often by only one male per spawning bout (Ref. 86779). |
| SoftBottom | -1 |
| HardBottom | -1 |
| Rocky | -1 |
| SeaGrassBeds | -1 |
| Entered | 2 |
| Dateentered | 1991-10-17T00:00:00.000Z |
| Modified | 2374 |
| Datemodified | 2014-02-06T00:00:00.000Z |
| SpecCode | 69 |
Taxize
We can also get it’s taxonomy rather easily. The package taxize allows you to extract and validate, among other things, the taxonomy of millions of species by accessing an important number of online databases accessible through their Application program interface (API).
- Encyclopedia of Life (EOL) eol
- axonomic Name Resolution Service tnrs
- Integrated Taxonomic Information Service (ITIS) itis
- Global Names Resolver (from EOL/GBIF) gnr
- Global Names Index (from EOL/GBIF)
- IUCN Red List
- Tropicos (from Missouri Botanical Garden)
- Theplantlist.org
- Catalogue of Life
- National Center for Biotechnology Information
- CANADENSYS Vascan name search API
- International Plant Names Index (IPNI)
- World Register of Marine Species (WoRMS)
- Barcode of Life Data Systems (BOLD)
- Pan-European Species directories Infrastructure (PESI)
- Mycobank
- National Biodiversity Network (UK)
- Index Fungorum
- EU BON
- Index of Names (ION)
- Open Tree of Life (TOL)
- World Register of Marine Species (WoRMS)
- NatureServe
# Start by exporting the taxonomy of our species
taxo <- taxize::classification(sp, db = 'worms', verbose = FALSE)
taxo[[1]]## # A tibble: 10 × 3
## name rank id
## <chr> <chr> <int>
## 1 Animalia Kingdom 2
## 2 Chordata Phylum 1821
## 3 Vertebrata Subphylum 146419
## 4 Gnathostomata Superclass 1828
## 5 Pisces Superclass 11676
## 6 Actinopterygii Class 10194
## 7 Gadiformes Order 10313
## 8 Gadidae Family 125469
## 9 Gadus Genus 125732
## 10 Gadus morhua Species 126436 # We can also extract the common or scientific names using sci2comm() & comm2sci(), respectively.
taxize::sci2comm(sp, db = 'itis')##
## Retrieving data for taxon 'Gadus morhua'## $`Gadus morhua`
## [1] "morue de l'Atlantique" "bacalao del Atl<e1>ntico"
## [3] "cod" "rock cod"
## [5] "morue franche" "Atlantic cod" # Or find out whether there are other names under which the species is known
taxize::synonyms(sp, db = 'itis')##
## Retrieving data for taxon 'Gadus morhua'## $`Gadus morhua`
## sub_tsn acc_tsn message
## 1 164712 164712 no syns found
##
## attr(,"class")
## [1] "synonyms"
## attr(,"db")
## [1] "itis" # Another really interesting feature is to extract all known species at a given taxonomic scale.
# Let's take the genus level and see the first 20 species that are part of that genus on the itis database
taxize::children(strsplit(sp, ' ')[[1]][1], db = 'itis')[[1]]##
## Retrieving data for taxon 'Gadus'## # A tibble: 4 × 5
## parentname parenttsn rankname taxonname tsn
## * <chr> <chr> <chr> <chr> <chr>
## 1 Gadus 164710 Species Gadus macrocephalus 164711
## 2 Gadus 164710 Species Gadus morhua 164712
## 3 Gadus 164710 Species Gadus ogac 164717
## 4 Gadus 164710 Species Gadus chalcogrammus 934083GloBI
How about extracting all known preys and predators of the species of interest? The Global Biotic Interactions web platform contains thousands of empirical binary interactions for multiple types of interactions, all over the world, and is easily accessible using the package rglobi.
# There are multiple types of interactions available on GloBI
rglobi::get_interaction_types()[,1:3]## interaction source target
## 1 eats consumer food
## 2 eatenBy food consumer
## 3 preysOn predator prey
## 4 preyedUponBy prey predator
## 5 kills killer victim
## 6 killedBy victim killer
## 7 parasiteOf parasite host
## 8 hasParasite host parasite
## 9 hostOf host symbiont
## 10 hasHost symbiont host
## 11 pollinates pollinator plant
## 12 pollinatedBy plant pollinator
## 13 pathogenOf pathogen host
## 14 hasPathogen host pathogen
## 15 vectorOf vector pathogen
## 16 hasVector pathogen vector
## 17 dispersalVectorOf vector seed
## 18 hasDispersalVector seed vector
## 19 symbiontOf source target
## 20 flowersVisitedBy plant visitor
## 21 visitsFlowersOf visitor plant
## 22 interactsWith source target # For now let's focus on predator-prey interactions
prey <- rglobi::get_prey_of(sp)$target_taxon_name
pred <- rglobi::get_predators_of(sp)$target_taxon_name
prey## [1] "Neocalanus tonsus" "Pseudocalanus elongatus"
## [3] "Arctica islandica" "Gastrosaccus spinifer"
## [5] "Diastylis rathkei" "Buccinum undatum"
## [7] "Corystes cassivelanus" "Eledone cirrhosa"
## [9] "Gonatus fabricii" "Bathypolypus arcticus"
## [11] "Rossia moelleri" "Bathypolypus bairdii"
## [13] "Todarodes sagittatus" "Cancer pagurus"
## [15] "Rossia macrosoma" "Pandalus borealis"
## [17] "Pandalus montagui" "Lithodes maja"
## [19] "Hyas coarctatus" "Crangon allmanni"
## [21] "Galathea strigosa" "Ophiopholis aculeata"
## [23] "Natatolana borealis" "Atelecyclus rotundatus"
## [25] "Crangon crangon" "Carcinus maenas"
## [27] "Mya arenaria" "Pagurus bernhardus"
## [29] "Ophiothrix fragilis" "Psammechinus miliaris"
## [31] "Trisopterus luscus" "Callionymus lyra"
## [33] "Actinopterygii" "Gadiformes"
## [35] "Gobiidae" "Eutrigla gurnardus"
## [37] "Trachurus trachurus" "Hippoglossoides platessoides"
## [39] "Scomber scombrus" "Lepidorhombus whiffiagonis"
## [41] "Trisopterus esmarkii" "Sardina pilchardus"
## [43] "Trisopterus minutus" "Enchelyopus cimbrius"
## [45] "Buglossidium luteum" "Microchirus variegatus"
## [47] "Merlangius merlangus" "Actinopterygii"
## [49] "Ammodytes tobianus" "Clupea harengus"
## [51] "Mallotus villosus" "Lethenteron camtschaticum"
## [53] "Lumpenus lampretaeformis" "Sprattus sprattus"
## [55] "Pomatoschistus minutus" "Zoarces viviparus"
## [57] "Gadus morhua" "Ammodytes marinus"
## [59] "Buenia jeffreysii" "Lophius piscatorius"
## [61] "Sebastes viviparus" "Microstomus kitt"
## [63] "Agonus cataphractus" "Solea solea"
## [65] "Merluccius merluccius" "Micromesistius poutassou"
## [67] "Salmo salar" "Myxine glutinosa"
## [69] "Glyptocephalus cynoglossus" "Scophthalmus rhombus"
## [71] "Pleuronectes platessa" "Pollachius virens"
## [73] "Scophthalmus maximus" "Zeugopterus punctatus"
## [75] "Melanogrammus aeglefinus" "Phrynorhombus norvegicus"
## [77] "Squalus acanthias" "Merluccius bilinearis"
## [79] "Zoarces americanus" "Pseudopleuronectes americanus"
## [81] "Scophthalmus aquosus" "Ammodytes dubius"
## [83] "Limanda limanda" "Myoxocephalus scorpius"
## [85] "Scyliorhinus canicula" pred## [1] "Larus" "Thalasseus sandvicensis"
## [3] "Phalacrocorax carbo" "Phoca vitulina"
## [5] "Sprattus sprattus" "Clupea harengus"
## [7] "Scomber scombrus" "no name"
## [9] "Gadus morhua" "Thunnus thynnus"
## [11] "Petromyzon marinus" "Pollachius virens"
## [13] "Molva molva" "Squalus acanthias"
## [15] "Hippoglossoides platessoides" "Pomatomus saltatrix"
## [17] "Myxine glutinosa" "Lophius piscatorius"
## [19] "Eutrigla gurnardus" "Hippoglossus hippoglossus"
## [21] "Amblyraja radiata" "Anarhichas lupus"
## [23] "Reinhardtius hippoglossoides" "Sebastes"
## [25] "Sebastes mentella" "Xiphias gladius"
## [27] "Merlangius merlangus"Make your search spatially explicit
Get data from OBIS
OBIS is the Ocean Biogeographic Information System and their vision is: “To be the most comprehensive gateway to the world’s ocean biodiversity and biogeographic data and information required to address pressing coastal and world ocean concerns.”
We can get access to their HUGE database through the excellent robis package by ropensci (unsurprisingly, they also have other great open science R packages, from which this whole section is highly inspired)
sf::st_as_text(bb$geometry)
OBIS <- robis::occurrence(scientificname = sp, geometry=st_as_text(bb$geometry), startdate = as.Date("2010-01-01"), enddate = as.Date("2017-01-01"), fields = c("species", "yearcollected","decimalLongitude", "decimalLatitude"))
# remove duplicates
OBIS <- unique(OBIS)
write.csv(OBIS,'obis.csv',row.names = FALSE)Visualize data
We can easily visualize the data…
OBIS <- read.csv("obis.csv")
plot(OBIS[,c("decimalLongitude", "decimalLatitude")])
… but the aesthetics could be much better!
Plot your study site
Define site map limits
Let’s start by defining a bounding box for your study site:
#create a matrix
bbmat <- cbind(
c(lonmin,lonmax,lonmax,lonmin,lonmin),
c(latmin,latmin,latmax,latmax,latmin)
)
# make the matrix a 'simple features' polygon
bbsf <- st_polygon(list(bbmat))
# and let's give it information about the projection:
bbsfc <- st_sfc(bbsf,crs="+proj=longlat +datum=WGS84")
# finally, let's make it a simple features data.frame
bb <- st_sf(name="Study Site",geometry=bbsfc)
str(bb)## Classes 'sf' and 'data.frame': 1 obs. of 2 variables:
## $ name : chr "Study Site"
## $ geometry:sfc_POLYGON of length 1; first list element: List of 1
## ..$ : num [1:5, 1:2] -71.1 -55.7 -55.7 -71.1 -71.1 ...
## ..- attr(*, "class")= chr "XY" "POLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
## ..- attr(*, "names")= chr "name"# plot it!
ggplot(bb) + geom_sf()
Add a basemap!
There are many ways to add a basemap in R
marmap package
bathymetry <- getNOAA.bathy(lonmin,lonmax,latmin,latmax,resolution=1,keep=TRUE)## File already exists ; loading 'marmap_coord_-71.06333;45.52399;-55.73636;52.01312_res_1.csv'plot.bathy(bathymetry,image=T,drawlabels=TRUE)
blues <- colorRampPalette(c("darkblue", "lightblue1"))
greys <- colorRampPalette(c(grey(0.4),grey(0.99)))
plot.bathy(bathymetry,
image = TRUE,
land = TRUE,
n=0,
bpal = list(c(0, max(bathymetry), greys(100)),
c(min(bathymetry), 0, blues(100))))
# this is a 'bathy' object and you could plot it as is with the tools provided in marmap, but to keep things 'simple' I'm going to convert to sf. There is no direct method yet, so... don't judge me!
# convert from bathy to raster to 'sp' polygon to simple features
# bathypoly <- marmap::as.raster(bathymetry) %>%
# rasterToPolygons() %>% st_as_sf
bathyras <- fortify.bathy(bathymetry)
bathyras$z[bathyras$z>0] <- NA
x=0.2
ggplot() +
geom_raster(data=bathyras,aes(x=x,y=y,fill=z))+
geom_sf(data=bb,fill="transparent",colour='yellow',size=2)+
# geom_sf(data=Canada,fill='grey40')+
coord_sf(xlim = c(lonmin-x,lonmax+x),
ylim = c(latmin-x,latmax+x),
expand = F)
Overlay your species occurrences over the base map
You can now overlay your species occurrences on your base map!
OBIS <- st_as_sf(OBIS,
coords = c("decimalLongitude", "decimalLatitude"),
crs="+proj=longlat +datum=WGS84",
remove=FALSE) %>%
filter(!is.na(species))
ggplot(OBIS) +
geom_raster(data=bathyras,aes(x=x,y=y,fill=z))+
geom_sf(data=bb,fill="transparent",colour='yellow',size=2)+
# geom_sf(data=Canada,fill='grey40')+
geom_point(data=OBIS,aes(x=decimalLongitude,y=decimalLatitude))+
facet_wrap(~species)+
coord_sf(xlim = c(lonmin-x,lonmax+x),
ylim = c(latmin-x,latmax+x),
expand = F)
Environmental variables
There are also multiple resources to access environmental data. As seen before, we can easily access bathymetry data using the marmap package. Other environmental data can also be accessed using other packages like raster for terrestrial climatic data and rnoaa
Extract sea surface temperature (SST)
# Bear with us, for simplicity's sake we will only take one year of data for one month, knowing very well that this does not make sense ecologically!
# Eventually use apply() to load multiple datasets and create rasterStack
# An exemple with sea ice
# urls <- seaiceeurls(mo='Apr', pole='N')[1:10]
# out <- lapply(urls, seaice)
# names(out) <- seq(1979,1988,1)
# df <- ldply(out)
# library('ggplot2')
# ggplot(df, aes(long, lat, group=group)) +
# geom_polygon(fill="steelblue") +
# theme_ice() +
# facet_wrap(~ .id)
sst.nc <- rnoaa::ersst(2010, 8, overwrite = TRUE)
sst <- raster(sst.nc$filename, varname = 'sst')
sst@data@values <- values(sst) # make sure values are in the right data slot
# never do that...
sst@extent@xmin <- sst@extent@xmin - 361
sst@extent@xmax <- sst@extent@xmax - 361
# Transform rasterLayer as a data.frame
sst.spdf <- as(sst, "SpatialPixelsDataFrame")
sst.df <- as.data.frame(sst.spdf)
colnames(sst.df) <- c('sst','x','y')
# Plot newly acquired data!
ggplot(sst.df, aes(x=x, y=y)) +
geom_tile(aes(x=x,y=y,fill=sst)) +
coord_equal()
Species distribution models
Now what we can do once we have species occurrence data and environmental variables are species distribution models (SDMs). There are multiple ways and packages to perform SDMs and this vignette provides an extensive overview of the different methods available. For this workshop, we will be using the package biomod2, which implements most of the methods presented in the vignette. There is also a thorough explanation of our to use biomod2 in another very interesting vignette
Formating data
The first thing to do is to format our data to be able to work with biomod2 functions, which work with BIOMOD.formated.data.PA class objects. This entails stacking all of our environmental data together to perform analyses.
# Create raster stack with environmental variables
# raster have to be of similar extent and resolution
# Load bathy data from marmap as raster
bathypoly <- marmap::as.raster(bathymetry)
bathypoly@data@values[bathypoly@data@values > 0] <- NA
# Resample sst to make it same extent as bathypoly
# create empty raster with appropriate resolution and extent
sst2 <- raster(resolution = res(bathypoly), ext = extent(bathypoly), crs=bathypoly@crs)
# Beware, this may not be good practice because we are changing the resolution of the data to a finer resolution!
sst2 <- resample(sst, sst2, method='bilinear') # extract data from sst
# Create raster stack
envCov <- raster::stack(bathypoly, sst2)
plot(envCov)
You can then create your biomod2 data class.
OBIS <- read.csv("obis.csv")
SDMdata <- biomod2::BIOMOD_FormatingData(resp.var = rep(1,nrow(OBIS)),
expl.var = envCov,
resp.xy = OBIS[, c("decimalLongitude","decimalLatitude")],
resp.name = sp,
PA.nb.rep = 2)##
## -=-=-=-=-=-=-=-=-=-=-=-= Gadus morhua Data Formating -=-=-=-=-=-=-=-=-=-=-=-=
##
## Response variable name was converted into Gadus.morhua
## ! No data has been set aside for modeling evaluation
## > Pseudo Absences Selection checkings...
## > random pseudo absences selection
## > Pseudo absences are selected in explanatory variables
## ! Some NAs have been automaticly removed from your data
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= SDMdata##
## -=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=
##
## sp.name = Gadus.morhua
##
## 1187 presences, 0 true absences and 1908 undifined points in dataset
##
##
## 2 explanatory variables
##
## layer Extended.reconstructed.sea.surface.temperature
## Min. :-525.0 Min. :11.33
## 1st Qu.:-220.0 1st Qu.:15.72
## Median : -93.0 Median :16.70
## Mean :-144.6 Mean :16.39
## 3rd Qu.: -53.0 3rd Qu.:17.18
## Max. : 0.0 Max. :18.15
##
##
## 2 Pseudo Absences dataset available ( PA1 PA2 ) with 967
## absences in each (true abs + pseudo abs)
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= plot(SDMdata)
You can then use this object with the biomod2 functions are generate a distribution model for your species of interest.
# Basic options for modelling
SDMoption <- biomod2::BIOMOD_ModelingOptions()
# SDM model
SDM <- biomod2::BIOMOD_Modeling(data = SDMdata,
models = "MAXENT.Tsuruoka",
model.options = SDMoption)##
##
## Loading required library...
##
## Checking Models arguments...## Warning in .Models.check.args(data, models, models.options, NbRunEval,
## DataSplit, : Models will run with 'defaults' parameters##
## Creating suitable Workdir...
##
## ! Weights where automaticly defined for Gadus.morhua_PA1 to rise a 0.5 prevalence !
## ! Weights where automaticly defined for Gadus.morhua_PA2 to rise a 0.5 prevalence !
##
##
## -=-=-=-=-=-=-=-=-=-=-= Gadus.morhua Modeling Summary -=-=-=-=-=-=-=-=-=-=-=
##
## 2 environmental variables ( layer Extended.reconstructed.sea.surface.temperature )
## Number of evaluation repetitions : 1
## Models selected : MAXENT.Tsuruoka
##
## Total number of model runs : 2
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
##
## -=-=-=- Run : Gadus.morhua_PA1
##
##
## -=-=-=--=-=-=- Gadus.morhua_PA1_Full # Print model summary
SDM##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-= BIOMOD.models.out -=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## Modeling id : 1493503397
##
## Species modeled : Gadus.morhua
##
## Considered variables : layer
## Extended.reconstructed.sea.surface.temperature
##
##
## Computed Models : Gadus.morhua_PA1_Full_MAXENT.Tsuruoka
## Gadus.morhua_PA2_Full_MAXENT.Tsuruoka
##
##
## Failed Models : none
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= # print the dimnames of this object to understand the structure of the returned objects
dimnames(SDM)## NULL # get all models evaluation.
# These give you an evaluation of the performance of your model.
biomod2::get_evaluations(SDM)## , , MAXENT.Tsuruoka, Full, PA1
##
## Testing.data Cutoff Sensitivity Specificity
## KAPPA 0.203 460.0 88.711 30.196
## TSS 0.193 460.0 88.711 30.196
## ROC 0.575 492.5 83.825 35.574
##
## , , MAXENT.Tsuruoka, Full, PA2
##
## Testing.data Cutoff Sensitivity Specificity
## KAPPA 0.190 455.0 90.480 27.311
## TSS 0.179 498.0 83.909 34.454
## ROC 0.567 497.5 83.909 34.454 # Finally, once the models are calibrated and evaluated, you can project your species spatial distribution...
SDMproj <- biomod2::BIOMOD_Projection(modeling.output = SDM,
new.env = envCov,
proj.name = 'Awesome!',
selected.models = 'all',
binary.meth = 'TSS',
compress = 'xz',
clamping.mask = F,
output.format = '.grd')##
## -=-=-=-=-=-=-=-=-=-=-=-=-= Do Models Projections -=-=-=-=-=-=-=-=-=-=-=-=-=
##
## > Building clamping mask
##
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer') # ... and visualize them!
plot(SDMproj)
Riveting! We starting with a basic idea and ended up with a species description, an idea of the distribution of certain environmental parameters in our area of interest (could be better, we know!) and the potential spatial distribution of our species!
…but should we stop there?
The reality is, there is so much more going on in our area of interest than the spatial distribution of a single species. Let’s see what other insights we an get from available open source resources.
Community level analyses
One of the usual challenges in ecology is to perform analyses at the community scale often due to the difficulty of sampling multiple species in a single system. The reality is that our species of interest is very likely to interact and be impacted by other species.
Extract all occurrence known occurrence data in area of interest
OBIS once again gives us the opportunity to access multiple species occurrences. Let’s go back to OBIS and extract every available data in our system.
# Load all known species occurrences from our area of interest
multiOBIS <- robis::occurrence(geometry=st_as_text(bb$geometry), startdate = as.Date("2010-01-01"), enddate = as.Date("2017-01-01"), fields = c("species", "yearcollected","decimalLongitude", "decimalLatitude"))
multiOBIS <- multiOBIS[!is.na(multiOBIS[,'species']), ] # remove species == NA
# Remove duplicates
multiOBIS <- unique(multiOBIS)
# Select species for which there is at least n observations
n = 50
spOK <- names(which(table(multiOBIS[,'species']) >= n)) # id of species with at least n observations
multiOBIS <- multiOBIS[multiOBIS[, 'species'] %in% spOK, ] # Subset of multiOBIS data with species with at least n observations
multiOBIS[,'species'] <- as.factor(multiOBIS[,'species'])
multiOBIS[,'yearcollected'] <- as.factor(multiOBIS[,'yearcollected'])
write.csv(multiOBIS,'multiOBIS.csv',row.names = FALSE) # Read occurrence data file
multiOBIS <- read.csv("multiOBIS.csv")
# Structure of object
str(multiOBIS)## 'data.frame': 37291 obs. of 4 variables:
## $ species : chr "Themisto compressa" "Themisto compressa" "Themisto compressa" "Themisto compressa" ...
## $ yearcollected : int 2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
## $ decimalLongitude: num -68.6 -68.2 -64 -63.7 -60.2 ...
## $ decimalLatitude : num 48.7 48.8 47.7 47.7 47 ... # Visualize species list
spList <- unique(multiOBIS[,'species'])
length(spList)## [1] 167 spList## [1] "Themisto compressa"
## [2] "Themisto abyssorum"
## [3] "Temora longicornis"
## [4] "Scina borealis"
## [5] "Scolecithricella minor"
## [6] "Parasagitta elegans"
## [7] "Paraeuchaeta norvegica"
## [8] "Triconia borealis"
## [9] "Oithona similis"
## [10] "Oithona atlantica"
## [11] "Metridia lucens"
## [12] "Metridia longa"
## [13] "Eukrohnia hamata"
## [14] "Dimophyes arctica"
## [15] "Calanus hyperboreus"
## [16] "Calanus glacialis"
## [17] "Calanus finmarchicus"
## [18] "Boreomysis arctica"
## [19] "Aglantha digitale"
## [20] "Acartia (Acartiura) longiremis"
## [21] "Reinhardtius hippoglossoides"
## [22] "Pandalus borealis"
## [23] "Myoxocephalus scorpius"
## [24] "Merluccius bilinearis"
## [25] "Hippoglossus hippoglossus"
## [26] "Glyptocephalus cynoglossus"
## [27] "Enchelyopus cimbrius"
## [28] "Chionoecetes opilio"
## [29] "Menidia menidia"
## [30] "Gadus morhua"
## [31] "Fundulus heteroclitus"
## [32] "Hyas araneus"
## [33] "Skeletonema costatum"
## [34] "Pseudo-nitzschia delicatissima"
## [35] "Homarus americanus"
## [36] "Carcinus maenas"
## [37] "Gasterosteus aculeatus"
## [38] "Apeltes quadracus"
## [39] "Gasterosteus wheatlandi"
## [40] "Palaemonetes vulgaris"
## [41] "Pseudopleuronectes americanus"
## [42] "Osmerus mordax"
## [43] "Fundulus diaphanus"
## [44] "Pungitius pungitius"
## [45] "Pleuronectes putnami"
## [46] "Cancer plebejus"
## [47] "Microgadus tomcod"
## [48] "Syngnathus fuscus"
## [49] "Tautogolabrus adspersus"
## [50] "Myoxocephalus aenaeus"
## [51] "Scophthalmus aquosus"
## [52] "Morone saxatilis"
## [53] "Urophycis tenuis"
## [54] "Alosa pseudoharengus"
## [55] "Brisaster fragilis"
## [56] "Ctenodiscus crispatus"
## [57] "Mallotus villosus"
## [58] "Triglops murrayi"
## [59] "Hyas coarctatus"
## [60] "Myoxocephalus octodecemspinosus"
## [61] "Leptoclinus maculatus"
## [62] "Illex illecebrosus"
## [63] "Aspidophoroides monopterygius"
## [64] "Liparis gibbus"
## [65] "Icelus spatula"
## [66] "Artediellus atlanticus"
## [67] "Hippoglossoides platessoides"
## [68] "Lycodes lavalaei"
## [69] "Limanda ferruginea"
## [70] "Gymnocanthus tricuspis"
## [71] "Artediellus uncinatus"
## [72] "Clupea harengus"
## [73] "Lithodes maja"
## [74] "Aspidophoroides olrikii"
## [75] "Amblyraja radiata"
## [76] "Ammodytes dubius"
## [77] "Eumesogrammus praecisus"
## [78] "Chlamys islandica"
## [79] "Lycodes vahlii"
## [80] "Hemitripterus americanus"
## [81] "Malacoraja senta"
## [82] "Arctozenus risso"
## [83] "Lumpenus lampretaeformis"
## [84] "Leptagonus decagonus"
## [85] "Eumicrotremus spinosus"
## [86] "Myxine glutinosa"
## [87] "Gadus ogac"
## [88] "Phycis chesteri"
## [89] "Gymnelus viridis"
## [90] "Scomber scombrus"
## [91] "Cryptacanthodes maculatus"
## [92] "Nezumia bairdii"
## [93] "Melanostigma atlanticum"
## [94] "Cyclopterus lumpus"
## [95] "Lophius americanus"
## [96] "Anarhichas lupus"
## [97] "Centroscyllium fabricii"
## [98] "Careproctus reinhardti"
## [99] "Lycenchelys verrillii"
## [100] "Phocoena phocoena"
## [101] "Boreogadus saida"
## [102] "Argentina silus"
## [103] "Pandalus montagui"
## [104] "Hippasteria phrygiana"
## [105] "Strongylocentrotus droebachiensis"
## [106] "Psilaster andromeda"
## [107] "Solaster endeca"
## [108] "Cucumaria frondosa"
## [109] "Crossaster papposus"
## [110] "Lebbeus polaris"
## [111] "Sebastes mentella"
## [112] "Pontophilus norvegicus"
## [113] "Oceanites oceanicus"
## [114] "Morus bassanus"
## [115] "Uria aalge"
## [116] "Uria lomvia"
## [117] "Fulmarus glacialis"
## [118] "Rissa tridactyla"
## [119] "Larus marinus"
## [120] "Alle alle"
## [121] "Larus argentatus"
## [122] "Alca torda"
## [123] "Fratercula arctica"
## [124] "Oceanodroma leucorhoa"
## [125] "Puffinus griseus"
## [126] "Puffinus gravis"
## [127] "Hemiselmis virescens"
## [128] "Quadricilia rotundata"
## [129] "Amphidinium sphenoides"
## [130] "Ceratoneis closterium"
## [131] "Lennoxia faveolata"
## [132] "Chaetoceros convolutus"
## [133] "Pyramimonas plurioculata"
## [134] "Telonema subtile"
## [135] "Plagioselmis prolonga"
## [136] "Leucocryptos marina"
## [137] "Pseudopedinella pyriformis"
## [138] "Heterocapsa rotundata"
## [139] "Gyrodinium fusus"
## [140] "Leptocylindrus minimus"
## [141] "Gyrodinium flagellare"
## [142] "Protoperidinium bipes"
## [143] "Lohmaniella oviformis"
## [144] "Gymnodinium verruculosum"
## [145] "Gyrodinium formosum"
## [146] "Gymnodinium galeatum"
## [147] "Torodinium robustum"
## [148] "Teleaulax amphioxeia"
## [149] "Chaetoceros similis"
## [150] "Monosiga marina"
## [151] "Rhodomonas marina"
## [152] "Peridiniella danica"
## [153] "Amphidinium kesslitzii"
## [154] "Gyrodinium spirale"
## [155] "Prorocentrum cordatum"
## [156] "Mesodinium rubrum"
## [157] "Katodinium glaucum"
## [158] "Chaetoceros debilis"
## [159] "Attheya septentrionalis"
## [160] "Pseudoscourfieldia marina"
## [161] "Gymnodinium elongatum"
## [162] "Thalassiosira nordenskioeldii"
## [163] "Nematopsides vigilans"
## [164] "Amphidinium crassum"
## [165] "Apedinella radians"
## [166] "Eutreptiella gymnastica"
## [167] "Thalassiosira pacifica" # Plot
multiOBIS <- sf::st_as_sf(multiOBIS,
coords = c("decimalLongitude", "decimalLatitude"),
crs="+proj=longlat +datum=WGS84",
remove=FALSE) %>%
filter(!is.na(species))
ggplot(multiOBIS) +
geom_raster(data=bathyras,aes(x=x,y=y,fill=z))+
geom_sf(data=bb,fill="transparent",colour='yellow',size=2)+
geom_point(data=multiOBIS,aes(x=decimalLongitude,y=decimalLatitude))+
coord_sf(xlim = c(lonmin-x,lonmax+x),
ylim = c(latmin-x,latmax+x),
expand = F)
Trophic network
nsp <- length(spList)
foodWeb <- matrix(ncol = nsp, nrow = nsp, data = 0, dimnames = list(spList, spList))
pb <- txtProgressBar(min = 0,max = nsp, style = 3)
for(i in 1:nsp) {
for(j in 1:nsp) {
inter <- rglobi::get_interactions_by_taxa(sourcetaxon = spList[i], targettaxon = spList[j], interactiontype = 'eats')
if(nrow(inter) > 0) foodWeb[j,i] <- 1
}
setTxtProgressBar(pb, i)
}
close(pb)
saveRDS(foodWeb, file = './foodWeb.rds')We can now visualize the network using the package iGraph
nsp <- length(spList)
foodWeb <- readRDS('./foodWeb.rds')
netw <- igraph::graph_from_adjacency_matrix(t(foodWeb))
vec_col <- spList
vec_col[vec_col != sp] <- "#11bdec"
vec_col[vec_col == sp] <- "#cc154b"
plot(netw,
vertex.color = vec_col,
vertex.frame.color = "transparent",
vertex.label.color = "black",
edge.arrow.size = .5
)
… Not very nice visually when there are too many species. We can however do this better (and more fun!) using the package networkD3
netwD3 <- networkD3::igraph_to_networkD3(netw, group = vec_col, what = 'both')
nodeSize <- spList
nodeSize[nodeSize != sp] <- 2
nodeSize[nodeSize == sp] <- 4
netwD3$nodes <- cbind(netwD3$nodes, nodeSize)
networkD3::forceNetwork(Links = netwD3$links,
Nodes = netwD3$nodes,
Source = 'source',
Target = 'target',
NodeID = 'name',
Group = 'group',
zoom = TRUE,
linkDistance = 50,
fontSize = 12,
opacity = 0.9,
charge = -20,
# colourScale = networkD3::JS('d3.scale.ordinal().range(vec_col)'),
Nodesize = 'nodeSize')Neat huh?!?
Co-distribution
We could also evaluate how species are co-distributed using an analysis analog to SDMs, but at the communicty scale and using a hierarchical Bayesian approach, Hierarchical Modeling of Species Communities (HMSC; ref) currently being implemented in R in the package HMSC (ref). As with SDM using the biomod2 package, the first step is to format the data into a HMSC class object.
Formatting the data
Create regular grid over study area
In order to create common observations among individual observations, we will aggregate all data within regular grid cells, so that cells will become planning units on which we base co-distribution of species. We will also create a rather coarse sized grid in order to make processing speed more efficient. Of course there would be more sophisticated ways to do this, but for simplicity’s sake we will perform a simple intersect between the grid and occurrence data.
# First create a function to generate hexagonal grids
hexGrid <- function(size, shp_over) {
library(rgeos)
library(sp)
shp_over2 <- gBuffer(shp_over, width = size) # buffer so that whole surface of shapefile is covered by a cell
spat_grid <- spsample(shp_over2,type="hexagonal",cellsize=size, offset=c(0,0)) # creating points for a hex grid
spat_grid <- HexPoints2SpatialPolygons(spat_grid) # creating polygons
proj4string(spat_grid) <- proj4string(shp_over) # set projection
spat_grid <- spat_grid[shp_over, ] # clipping grid to retain only cells intersecting shapefile
id <- paste('ID',seq(1:length(spat_grid)), sep = '')
row.names(spat_grid) <- id
spat_grid <- SpatialPolygonsDataFrame(Sr = spat_grid, data = data.frame(ID = id, row.names = id))
return(spat_grid)
}
# Load shapefile over which grid is made
egsl <- rgdal::readOGR(dsn = "./", layer = "egsl")
# Create grid
egsl_grid <- hexGrid(size = 30000, shp_over = egsl)
rgdal::writeOGR(egsl_grid, dsn = './', layer = 'egsl_grid', driver="ESRI Shapefile", overwrite_layer=TRUE)
# Plot grid
plot(egsl_grid) # Plot grid
egsl_grid <- rgdal::readOGR(dsn = "./", layer = "egsl_grid")## OGR data source with driver: ESRI Shapefile
## Source: "./", layer: "egsl_grid"
## with 419 features
## It has 1 fields proj <- "+proj=longlat +datum=WGS84"
egsl_grid <- sp::spTransform(egsl_grid, CRSobj = CRS(proj))
plot(egsl_grid)
Intersect occurrences and environmental data with grid
# reload multiOBIS data and transform into sp points
multiOBIS <- read.csv("multiOBIS.csv")
proj <- "+proj=longlat +datum=WGS84"
multiOBISpt <- sp::SpatialPointsDataFrame(coords = multiOBIS[, c("decimalLongitude", "decimalLatitude")], data = multiOBIS, proj4string = CRS(proj), match.ID = TRUE)
# Overlay occurrences on grid
occGrid <- sp::over(multiOBISpt, egsl_grid)
multiOBIS <- cbind(multiOBIS, occGrid)
multiOBIS_wide <- reshape2::dcast(multiOBIS, formula = ID ~ species, value.var = 'species', fun.aggregate = length)
for(i in 2:ncol(multiOBIS_wide)) {
for(j in 1:nrow(multiOBIS_wide)) {
if(multiOBIS_wide[j,i] > 1) multiOBIS_wide[j,i] <- 1
}
}
toAdd <- egsl_grid@data[,'ID'][!egsl_grid@data[,'ID'] %in% multiOBIS[,'ID']]
matSup <- matrix(nrow = length(toAdd), ncol = ncol(multiOBIS_wide), dimnames = list(c(), colnames(multiOBIS_wide)), data = 0)
matSup[,'ID'] <- toAdd
multiOBIS_wide <- rbind(multiOBIS_wide, matSup)
multiOBIS_wide <- multiOBIS_wide[-which(is.na(multiOBIS_wide[,'ID']) == T), ]
egsl_grid@data <- merge(egsl_grid@data, multiOBIS_wide, by = 'ID')
sp <- colnames(egsl_grid@data)[2:ncol(egsl_grid@data)]
spID <- seq(1:length(sp))
colnames(egsl_grid@data)[2:ncol(egsl_grid@data)] <- spID
# Extract environmental values from rasters
egsl_grid@data <- cbind(egsl_grid@data, raster::extract(bathypoly, egsl_grid, fun = mean, na.rm = T))
colnames(egsl_grid@data)[ncol(egsl_grid@data)] <- 'bathy'
egsl_grid@data <- cbind(egsl_grid@data, raster::extract(sst2, egsl_grid, fun = mean, na.rm = T))
colnames(egsl_grid@data)[ncol(egsl_grid@data)] <- 'sst'
# HMSC data
studyGrid <- egsl_grid@data
envCov <- c('bathy','sst')
# Remove environmental NAs
NAs <- studyGrid[, envCov] %>%
lapply(X = ., FUN = function(x) which(is.na(x))) %>%
unlist(.) %>%
unique(.)
studyGrid <- studyGrid[-NAs, ]Create HMSC class data
# ---------------------------------
# Y: sample units by species matrix
# ---------------------------------
# Extracting only presence/absence data
Station <- studyGrid[, 'ID']
Y <- studyGrid[, as.character(spID)]
rownames(Y) <- Station
for(i in 1:ncol(Y)) {
Y[,i] <- as.numeric(as.character(Y[,i]))
}
# -----------------------------------------
# Pi: sample units by random effects matrix
# -----------------------------------------
# Create a dataframe for random effects, columns have to be factors
Pi <- data.frame(sampling_unit = as.factor(Station))
# ----------------------------------------------------
# X: sampling units by environmental covariates matrix
# ----------------------------------------------------
# Create a matrix for the values of environmental covariates at each sampling unit location
X <- matrix(ncol = length(envCov),
nrow = nrow(studyGrid))
X <- studyGrid[, envCov]
rownames(X) <- Station
# ----------------------------
# as.HMSCdata for HMSC package
# ----------------------------
# Creating HMSC dataset for analyses
HMSCdata <- HMSC::as.HMSCdata(Y = Y, X = X, Random = Pi, interceptX = T, scaleX = T)## [1] "column names were added to 'Random'"
## [1] "'Y' was converted to a matrix" str(HMSCdata)## List of 3
## $ Y : num [1:369, 1:167] 0 0 0 0 1 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:369] "ID10" "ID100" "ID101" "ID103" ...
## .. ..$ : chr [1:167] "1" "2" "3" "4" ...
## $ X : num [1:369, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:369] "ID10" "ID100" "ID101" "ID103" ...
## .. ..$ : chr [1:3] "Intercept" "bathy" "sst"
## $ Random:'data.frame': 369 obs. of 1 variable:
## ..$ random1: Factor w/ 369 levels "ID10","ID100",..: 1 2 3 4 5 6 7 8 9 10 ...
## - attr(*, "class")= chr "HMSCdata"HMSC analysis
We are now ready to perform the analyses!
model <- HMSC::hmsc(HMSCdata,
family = "probit",
niter = 10000,
nburn = 100,
thin = 10)
saveRDS(model, file = './HMSCmodel.rds')HMSC analysis diagnostics
model <- readRDS('./HMSCmodel.rds')
# =========================================
# 3. Producing MCMC trace and density plots
# =========================================
mixingMeansParamX <- coda::as.mcmc(model, parameters = "meansParamX")
# Trace and density plots to visually diagnose mcmc chains
# Another way to check for convergence is to use diagnostic tests such as Geweke's convergence diagnostic (geweke.diag function in coda) and the Gelman and Rubin's convergence diagnostic (gelman.diag function in coda).
paramModel <- colnames(mixingMeansParamX)
nParam <- length(paramModel)
par(mfrow = c(nParam, 2), mar = rep(2, 4))
for(i in 1:ncol(mixingMeansParamX)) {
coda::traceplot(mixingMeansParamX[,i], col = "blue", main = paste('Trace of ', paramModel[i]))
coda::densplot(mixingMeansParamX[,i], col = "orange", main = paste('Density of ', paramModel[i]))
}
# ================================
# 4. Producing posterior summaries
# ================================
# Average
average <- apply(model$results$estimation$paramX, 1:2, mean)
# 95% confidence intervals
CI.025 <- apply(model$results$estimation$paramX, 1:2, quantile, probs = 0.025)
CI.975 <- apply(model$results$estimation$paramX, 1:2, quantile, probs = 0.975)
# Summary table
paramXCITable <- cbind(unlist(as.data.frame(average)),
unlist(as.data.frame(CI.025)),
unlist(as.data.frame(CI.975)))
colnames(paramXCITable) <- c("average", "lowerCI", "upperCI")
rownames(paramXCITable) <- paste(rep(colnames(average), each = nrow(average)), "_", rep(rownames(average), ncol(average)), sep="")
# Credible intervals
paramXCITable_Full <- paramXCITable
beg <- seq(1,nrow(paramXCITable_Full), by = nsp)
end <- seq(nsp,nrow(paramXCITable_Full), by = nsp)
sign <- abs(as.numeric(paramXCITable_Full[, 'lowerCI'] <= 0 & paramXCITable_Full[, 'upperCI'] >=0) - 3)
# Export figure
par(mfrow = c((length(beg)+1),1))
for(i in 1:length(beg)) {
paramXCITable <- paramXCITable_Full[beg[i]:end[i], ]
cols <- sign[beg[i]:end[i]]
par(mar=c(1,2,1,1))
plot(0, 0, xlim = c(1, nrow(paramXCITable)), ylim = round(range(paramXCITable)), type = "n", xlab = "", ylab = "", main=paste(colnames(mixingMeansParamX)[i]), xaxt="n", bty = 'n')
axis(1,1:nsp,las=2, labels = rep('',nsp))
abline(h = 0,col = 'grey')
arrows(x0 = 1:nrow(paramXCITable), x1 = 1:nrow(paramXCITable), y0 = paramXCITable[, 2], y1 = paramXCITable[, 3], code = 3, angle = 90, length = 0.05, col = cols)
points(1:nrow(paramXCITable), paramXCITable[,1], pch = 15, cex = 1, col = cols)
}
mtext(text = sp, side = 1, line = 1, outer = FALSE, at = 1:nsp, col = 1, las = 2, cex = 0.4)
# =======================
# 6. Association networks
# =======================
# Extract all estimated associatin matrix
assoMat <- HMSC::corRandomEff(model)
# Average
plotMean <- apply(assoMat[, , , 1], 1:2, mean)
colnames(plotMean) <- rownames(plotMean) <- sp
# Build matrix of colours for chordDiagram
plotDrawCol <- matrix(NA, nrow = nrow(plotMean), ncol = ncol(plotMean))
plotDrawCol[which(plotMean > 0.4, arr.ind=TRUE)]<-"red"
plotDrawCol[which(plotMean < -0.4, arr.ind=TRUE)]<-"blue"
# Build matrix of "significance" for corrplot
plotDraw <- plotDrawCol
plotDraw[which(!is.na(plotDraw), arr.ind = TRUE)] <- 0
plotDraw[which(is.na(plotDraw), arr.ind = TRUE)] <- 1
plotDraw <- matrix(as.numeric(plotDraw), nrow = nrow(plotMean), ncol = ncol(plotMean))
# plots
# Matrix plot
Colour <- colorRampPalette(c("red", "white", "blue"))(200)
corrplot::corrplot(plotMean, method = "color", col = Colour, type = "lower", diag = FALSE, p.mat = plotDraw, tl.srt = 65, tl.cex = 0.4, pch.cex = 0.3, tl.col = 'black')
# Chord diagram
circlize::chordDiagram(plotMean, symmetric = TRUE, annotationTrack = 'grid', grid.col = "grey", col = plotDrawCol, preAllocateTracks = 1)
circlize::circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
xlim = circlize::get.cell.meta.data("xlim")
ylim = circlize::get.cell.meta.data("ylim")
sector.name = circlize::get.cell.meta.data("sector.index")
circlize::circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5), cex = 0.4)
}, bg.border = NA)
# ===============================================
# 7. Computing the explanatory power of the model
# ===============================================
# Prevalence
prevSp <- colSums(HMSCdata$Y)
# Coefficient of multiple determination
R2 <- HMSC::Rsquared(model, averageSp = FALSE)
# Draw figure
plot(prevSp, R2, xlab = "Prevalence", ylab = expression(R^2), cex = 0.8, pch=19, las=1, cex.lab = 1, main = 'Explanatory power of the model')
abline(h = mean(R2, na.rm = T), col = "blue", lwd = 2)
# =============================================
# 8. Generating predictions for validation data
# =============================================
source('./crossValid.r')
modelAUC <- crossValidation(data = HMSCdata, nCV = 2, validPct = 0.2, as.strata = FALSE)## [1] "The priors for the latent variables should be OK for probit models but not necessarily for other models, be careful"
## iteration 2000
## iteration 4000
## iteration 6000
## iteration 8000
## [1] "The priors for the latent variables should be OK for probit models but not necessarily for other models, be careful"
## iteration 2000
## iteration 4000
## iteration 6000
## iteration 8000 meanAUC <- colMeans(modelAUC)
sdAUC <- apply(modelAUC, MARGIN = 2, FUN = sd)
plot(prevSp, meanAUC, ylim = c(0,1), xlab = "Prevalence", ylab = 'AUC', cex = 0.8, pch=19, las=1, cex.lab = 1, main = 'Monte Carlo cross-validation with AUC of ROC curves')
abline(h = 0.5, col = 'grey')
arrows(x0 = prevSp, x1 = prevSp, y0 = (meanAUC - sdAUC), y1 = (meanAUC + sdAUC), code = 3, angle = 90, length = 0.05)## Warning in arrows(x0 = prevSp, x1 = prevSp, y0 = (meanAUC - sdAUC), y1 =
## (meanAUC + : zero-length arrow is of indeterminate angle and so skipped points(prevSp, meanAUC, pch = 19, cex = 0.8)